An Assessment of the State-of-the-art in 
Multidisciplinary Aeromechanical Analyses 

Anubhav Datta 
ELORET Corporation 
Ames Research Center 
Moffett Field, CA 
datta@merlin.arc.nasa.gov 


Wayne Johnson 
Aeromechanics Branch 
NASA Ames Research Center 
Moffett Field, CA 
wayne.johnson@nasa.gov 


ABSTRACT 

This paper presents a survey of the current state-of-the-art in multidisciplinary aeromechanical analyses 
which integrate advanced Computational Structural Dynamics (CSD) and Computational Fluid Dynamics 
(CFD) methods. The application areas to be surveyed include fixed wing aircraft, turbomachinary, and 
rotary wing aircraft. The objective of the authors in the present paper — together with a companion paper 
on requirements — is to lay out a path for a High Performance Computing (HPC) based next generation 


comprehensive rotorcraft analysis. From this survey 
possible to identify the critical technology gaps that 


INTRODUCTION 

This paper presents a survey of computational aeroe- 
lasticity in the disciplines of fixed wing aircraft, turboma- 
chinery, and rotary wing aircraft. The work was under- 
taken by the U.S. Army Aeroflightdynamics Directorate 
under the High Performance Computing Institute for Ad- 
vanced Rotorcraft Modeling and Simulation (HI-ARMS) 
and NASA. 

The survey covers the emerging methods which in- 
tegrate Reynolds-averaged Navier-Stokes (RANS) CFD, 
Finite Element Method (FEM) based structural mechan- 
ics, and high-fidelity coupling procedures designed to sat- 
isfy the unique requirements of each discipline. In each 
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of the key technologies in other application areas it is 
stem from unique rotorcraft requirements. 

discipline, the key aeromechanical phenomena which de- 
termine the costs and risks associated with design, but 
remain beyond current prediction capabilities, are de- 
scribed. The current status of High Performance Com- 
puting (HPC) based high-fidelity studies on the key phe- 
nomena are reviewed, the recent advances summarized, 
and the unresolved challenges highlighted. 

The central theme of this paper is rotorcraft. The 
objective of the authors — together with a companion pa- 
per on requirements [1] — is to lay out a path for a High 
Performance Computing (HPC) based truly comprehen- 
sive next generation rotorcraft code; comprehensive in 
solutions (performance, loads, stability, vibration, han- 
dling qualities), comprehensive in applications (ground, 
hover, steady flight, unsteady maneuvers), and compre- 
hensive in scope (arbitrary geometries, innovative con- 
figurations). The intention of the present review is to 
identify the key technologies in other application areas 
that can be drawn upon to this end, and to identify the 
critical technology gaps that stem from unique rotary 
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wing requirements. 

The paper is divided into five sections. The first 
section presents a description of the state-of-the-art in 
high-fidelity fluid structure coupling. This is followed by 
three sections, one each on the status of computational 
aeroelasticity in fixed wing aircraft, turbomachinery, and 
rotary wing aircraft. Each section is subdivided into two 
parts; the first part is on structural mechanics, the sec- 
ond part is on coupled fluid-structure applications. The 
CFD methods applicable to each discipline are not re- 
viewed in this paper. The last section summarizes the 
different CFD/CSD coupling nomenclatures used in the 
three disciplines. 

FLUID-STRUCTURE COUPLING 

Definition of the problem 

The problem of fluid-structure coupling involves 
three issues: (1) temporal accuracy, (2) spatial accuracy, 
and (3) interface geometry representation. The purpose 
of the present discussion is to clarify their meaning, ex- 
plain why they were of little importance in rotorcraft so 
far, and highlight why they are important now and for 
the future. 

At the PDE level, the three issues are related to the 
following questions. The first issue, that of temporal ac- 
curacy, is related to the question whether two systems 
of coupled PDEs of different types can be solved sepa- 
rately, one after the other, one step at a time, exchang- 
ing solutions at each time step. This is the method of 
partitioned formulations that begins with an acceptance 
that fluid PDEs of convection-diffusion type, and struc- 
tural PDEs of elliptic type are best solved separately in 
their own domains using their own efficient solvers. The 
main task then is to devise a method of solution exchange 
which renders the process at least as time accurate as the 
worse of the two solvers. The second and the third issues, 
those of spacial accuracy and interface geometry repre- 
sentation, are related to errors introduced during solution 
exchange across domains. Depending on these errors the 
method of solution exchange must be re-constructed to 
achieve an intended temporal accuracy. Note that the 
second and the third issues are independent of the first, 
that is, they arise whether or not a partitioned procedure 
is adopted. 

Partitioned vs. fully coupled formulations 

The practical benefits of partitioned formulations 
are obvious — modularity of framework, and domain sep- 
aration for refinements. However, an appropriate method 
of solution exchange must be designed to ensure time ac- 
curacy. Fully coupled formulations, on the other hand, 
are strictly time accurate without this additional bur- 


den. However, such formulations require: (1) a com- 
mon time integrator, and (2) solution of an algebraic sys- 
tem, at each Newton-like step, that includes both fluid 
and structural stiffness. The second requirement is not 
easy to meet for the following reasons. Direct solution 
is not an option, since modern CFD (without the use 
of reduced order modeling) provides far too many de- 
grees of freedom (DOFs), typically 10-100M for rotor- 
craft. Iterative solution is the natural option in CFD. 
The temporal evolution of a convection-diffusion equa- 
tion is naturally analogous to the iterative convergence 
of an algebraic system. Moreover, they are easy to par- 
allelize. On the other hand, iterative solvers are tradi- 
tionally not preferred for structures. The convergence 
rates of iterative solvers depend on the condition num- 
ber of the system (for symmetric positive definite ma- 
trices the condition number is the ratio of the largest to 
the smallest eigenvalues, or natural frequencies squared). 
Typically, aerospace structures pose 4-th order elastic- 
ity problems involving bending-torsion-extension of thin 
beams, or bending-shear-extension of thin plates, and 
shells. Condition numbers for such structures range from 
10 8 to 10 9 (the number for a typical rotor beam model 
is around 10 9 ). In comparison, the maximum condition 
number for a 2-nd order fluid problem can be as high as 
10 6 . Pre-conditioners that are well suited for structures 
are not suitable for fluids. For example, the incomplete 
Cholesky pre-conditioner cannot be applied to fluids as 
the fluid system is not symmetric positive definite. Simi- 
larly, pre-conditioners well suited for fluids, like the block 
Jacobi (easily parallelizable), can only be used for struc- 
tures with special re-ordering and sparsity fill-in for rea- 
sonable, yet problem dependent, convergence. Moreover, 
such procedures demonstrate poor scalability. Thus, the 
challenge of a unified solution procedure in the presence 
of CFD is real [2]. 

The difficulty is bypassed for structures that are less 
ill-conditioned. One example is turbomachinery struc- 
tures. These are modeled using solid elements, gov- 
erned by 2-nd order elasticity. Gottfried and Fleeter 
[3] have used full coupling for turbomachinery aeroelas- 
ticity. Their analyses, TAM-ALE3D, a refined version 
of the original ALE3D code developed in the Lawrence 
Livermore National Laboratories, is a fully coupled finite 
element analysis that has been applied to turbomachin- 
ery flutter calculations. A second example is biological 
structures. Here, structures are either membrane-like, or 
thick, when shell-like. Bathe and Zhang [4] have applied 
full coupling in biomedical fluid flows, as well as general 
internal flows in mechanical components. Their analysis 
is used extensively for biomedical and mechanical appli- 
cations, as part of the ADINA commercial software code. 
An option for partitioned formulation is also provided 
(the authors denote them as direct and iterative proce- 
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dures). The foundations of all fully coupled approaches 
can be traced back to the Arbitrary Lagrangian Eulerian 
(ALE) formulation of conservations laws, first introduced 
by Hirt et al. [5] in 1974 for fluid flows at all speeds. It 
was applied by Donea [6] in 1983 for fluid-structure in- 
teraction problems. Bendikson [7, 8] first applied this 
technique in the early 1990s for flutter calculations and 
demonstrated minimal errors in energy transfer between 
fluids and structures. Fully coupled formulations are also 
termed monolithic formulations. 

In summary, for aerospace structures, partitioned 
formulations provide fundamental advantages over fully 
coupled ones, in addition to their obvious practical ben- 
efits. 

Current coupling practices in rotorcraft 

Classical comprehensive analyses, which use lifting- 
line or lifting-surface models, involve at the most 1-10K 
DOFs. The formulation is fully coupled, and is solved 
using direct methods. Thus the first issue of temporal 
accuracy is satisfied. To re-iterate, the need for a par- 
titioned formulation for fluids and structures is felt only 
under the following two circumstances. First, when mil- 
lions of fluid DOFs call for an iterative solver whereas the 
structural DOFs, being ill-conditioned, call for a direct 
solver, or at least a different pre-conditioner. Second, 
when these conflicting requirements of fluids and struc- 
tures are very easily met in separate domains. The sec- 
ond issue of spatial accuracy is easily satisfied by such for- 
mulations as the structural shape functions are available 
to the fluid terms. The third issue of interface geometry 
representation has also been simple, so far, in rotorcraft. 
This is because the beam theory, or any single compo- 
nent structural model, regardless of 1-D, 2-D, or 3-D, 
provides a simple yet rigorous definition of the surface. 
The complexity arises in multi-component structures like 
the fuselage. 

Dynamics researchers have carried out coupled 
rotor-fuselage analysis, but at a time when CFD air- 
loads were beyond state of the art. With the emergence 
of this capability, there is a requirement to address the 
issue of geometry representation. A detailed dynamic 
model is not necessarily the best suited for fluid-structure 
coupling. The surface geometry representation is more 
important than the internal load paths. Thus, the re- 
quirement is to have, at the least, a shell model that 
re-produces the low frequency modes (up to 40 Hz), and 
at best, a detailed model that includes the outer skin. 
Fixed-wing researchers have already accomplished these 
tasks, and there is a volume of literature that can be 
drawn upon. 

Rotorcraft CFD/CSD coupled analyses have used 
partitioned procedures. The focus so far has been on the 


rotor. For an isolated rotor, the wetted surface is rigor- 
ously defined. The shape functions are available. Thus 
the issues of spatial accuracy and surface representation 
are easily dealt with. The issue of temporal accuracy is 
enforced concurrently with Newton-Raphson type itera- 
tions for determining the trim angles. For a time march- 
ing solution, the procedure is conceptually simpler as the 
trim angles are left undetermined, but requires the en- 
forcement of temporal accuracy at every time step. The 
problem is then the same as that faced by fixed-wing re- 
searchers, where it has been of great importance due to 
the emphasis on flutter. Small errors in coupling result in 
erroneous energy exchange and affect flutter boundaries 
adversely. The emphasis in rotors has been on loads. 
Errors in fluid-structure coupling are less visible. The 
problem is further compounded in fixed-wing because of 
its 3-D structure and multiple sub-structures. Without a 
detailed 3-D structure, this difficulty has not been faced 
yet in rotorcraft. 

Temporal accuracy 

The temporal accuracy of partitioned formulations 
is the main focus of high-fidelity fluid-structure coupling. 
Physically, it means that the airloads and structural 
loads at a given time step are consistent with one an- 
other. Numerically, it means that the temporal error is 
driven down to the level of the worse solver. 

One effective approach is the use of sub-iterations 
between every consecutive pair of time steps. The 
method was first applied in aeroelastic computations by 
Strganac and Mook [9] , followed by Weeratunga and Pra- 
mono [10], and more recently by Melville and his co- 
researchers [11, 12]. The method is time intensive for 
RANS. Potential innovations may involve the coupling 
of structural dynamics with the fluid sub-iterations. 

Early work in fixed-wing fluid structure coupling in- 
volved sub-iteration free, straight-forward time integra- 
tion. See for example, Edwards [13], Bennet [14], and 
Gurus wamy et al. [15, 16]. In the serial approach, one 
solver waited for the completion of the other, before ex- 
changing solutions and advancing to the next time step. 
In the parallel approach, both solvers advanced simulta- 
neously, followed by solution exchange before advancing 
to the next step. The latter was proposed originally by 
Weeratunga and Pramono [10], and refined subsequently 
by Farhat and Lesoinne [17]. The original work involved 
accompanying sub-iterations due to the relative instabil- 
ity of the scheme. The subsequent refinement improved 
stability without sub- iterations, at the cost of exchang- 
ing solutions twice at each time interval. Since then, a 
significant amount of research has been conducted on de- 
vising sub-iteration free methods that retain, at the least, 
second-order accuracy. Farhat and his co-researchers [18] 
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have described the formal design of such methods. 

Sub-iteration free approaches are predictor-corrector 
schemes, tailored to the individual time integrators. A 
generalized predictor-corrector approach is illustrated 
below, adopted from the last reference. The mesh defor- 
mation is denoted by x, the fluid variables are denoted by 
q , and the structural deformations by u. The time step 
is denoted by n. The fluid-structure interface boundary 
is denoted by T. For example Mr denotes the structural 
deformations at the interface, a subset of u. The symbol 
< — denotes a time integrator and shows dependencies 
on the latest time steps. A single time step advancement 
is given as follows. 

1. Predict deflections for the next time step 


5. Update structural solver 

u n+1 «— u n ,Ff x (including uf 1 ) 

The predictor is on displacement (step 1). The corrector 
is on forcing (step 4). Methods of the above type are 
classified as loose coupling in the fixed wing community. 
If the above procedure is carried out multiple times at the 
same time step (i.e. sub-iterations), then the following 
can be enforced 

uf lP = uf 1 

which implies that the mesh deformation follows 

X ^ 1 = X™ + u r ( U £ +1 - w?) 


u 


n+ 1 

r 


p 


2. Update mesh points x. Mesh boundary points de- 
noted by xt ■ Mesh internal points denoted by xq. 
The boundary mesh points can be updated as 


„n+i 


x r 


,n+l 


For example, 


xf lP = xf + U r (u” +lP - u£ P ) 


The above equation simply represents the discretized 
deflection and velocity compatibilities on the sur- 
face. Without re-griding in time, the deflection com- 
patibility is given by x = Upu where Ur is simply 
the transformation that connects the CFD surface 
grid to the structural grid. For perfectly matched 
meshes Ur is equal to the identity matrix. The ve- 
locity compatibility is x = Upu. If the structure 
reaches out to the wetted surface, and the shape 
functions are known, Up = H, where H are the 
shape functions. With re-griding, the mesh connec- 
tivity changes with time, i.e. the transformation Ur 
is now a function of time, and the discretized veloc- 
ity compatibility is expressed as x = Upii where Up 
denotes a selected combination of Up over previous 
time steps, e.g. the mean of Up and Up -1 . 


The internal points are updated based on the bound- 
ary points 


n+l 


T »+ 1 
> x r 


3. Update fluid solver using updated mesh 
q n+l < q n ,x n+lP 


rather than 

a-p+ lP = xf + Up (uP +lP - up") 

and hence enforces the velocity compatibility x = Uup 
strictly. The forcing and structural response (step 4 and 
5) are also consistent, as xf 1 is replaced with xf 1 . 
The method with sub-iterations is classified as close cou- 
pling, tight coupling, or strong coupling in the fixed wing 
literature. The force predictor in step 4 can be chosen in 
various ways. The first form is also called serial stagger- 
ing. The second form is called parallel staggering, simply 
because the fluid and structural updates (step 3 and 5) 
can be performed independently of one another. 

The methods that enforce second-order accuracy, 
without sub- iterations, rely on the formal selection of: 
(1) the integrators, i.e. the arrows — ’, (2) the initial 
predictor, i.e. step 1, and (3) the form of the force cor- 
rector, i.e. step 4. Two illustrative examples are provide 
in reference [18] based on a second-order fluid integra- 
tor, and a second-order (Newmark type) structural in- 
tegrator. Unless formally selected, the sub-iteration free 
methods provide only first-order accuracy. 

Lastly, we note that for non-CFD Lagrangian aero- 
dynamic models, the time iteration procedure is simpler 
as there are no grid motions. The need to reconcile the 
Geometric Conservation Law and structural motion up- 
date does not arise. The structural forcing F s is related 
directly to the deformations u. The serial and parallel 
staggered schemes along with sub-iterations can be im- 
plemented as follows. 

1. Serial staggered with sub- iterations 
Predict u n+lP 


4. Calculate structural forcing 


pn+l 


= F s (q n+1 ,xf lP ) 

or F s (q n ,x p V.. or other forms. 


Then, perform successively 

Ff 1 < — F™, u n+lP ; in step 1 

u n+ 1 < — m",F" +1 ; in step 2 

Iterate until convergence, i.e. u n+1 ss u n+1 
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2. Parallel staggered with sub-iterations 
Predict u n+1 

Then, perform in a single step 

pn+ 1 , pra , ra+1 . 

1 S ± S 5 U 5 

u n+1 < — u n ,F£\ 

Exchange updates, i.e. update u n+lP with u n+1 in 
the airloads integrator, and update F" with F” +1 
in the structures integrator, and iterate until con- 
vergence. 

The parallel staggered method is advantageous when 
the fluid and the structure are run on separate processors. 
The parallel method has inferior stability and requires 
requires lower time steps. 

Spatial accuracy 

Assume that the structural model includes a rig- 
orous interface geometry representation. However, the 
CFD and the CSD discretizations will not, in general, 
match at the interface. The deformation and struc- 
tural forcing must be mapped correctly across the non- 
matching interface. Generic methods which are ‘exact’ 
regardless of the blade shape are critical for the evalu- 
ation of advanced geometry blades. Two such methods 
are formulated below. 

Deformation mapping deals with the correct evalua- 
tion of Ur hr step 2 (see previous section). Loads transfer 
deals with the correct evaluation of structural forcing F s 
in step 4. The latter implies an exact evaluation of the 
virtual work term. This is a necessary energy conserva- 
tion requirement. In addition, it is desirable, though not 
necessary, that the total integrated forcing be preserved 
during the mapping. 

Deformation mapping is straight-forward when the 
shape functions are available. Loads transfer can be ac- 
complished via: (1) integrated force coupling, and (2) di- 
rect traction coupling. In the integrated force coupling, 
the virtual work is calculated based on integrated fluid 
stresses (pressure and skin friction) over a surface patch. 
In surface traction coupling, the virtual work is calcu- 
lated based on the pressure and shear distributions di- 
rectly. The first method preserves total forces, but does 
not guarantee energy conservation for a finite mesh. The 
second method is strictly energy conservative, but does 
not guarantee preservation of total forces. Thus, they 
are complimentary to each other. However, as long as in- 
terpolations and integrations are performed consistently 
within each domain, both satisfy and conservation and 
preservation in the limit of mesh refinement. 

Both methodologies can be formulated for rotor 
blades. Simple illustrations are given below. 
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Figure 1: A concentrated force on a rotor blade ob- 
tained by integrating fluid stresses over a surface 
patch of area A A 

The integrated force coupling is formulated as fol- 
lows. Consider a concentrated force F acting at a point 
P on the deformed blade, Fig. 1. F is integrated trac- 
tion over an arbitrary patch of area A A. The virtual 
work done by the force is simply 

SW = F ■ Sr P (1) 

where fp is the position vector of the point P, and 5fp is 
its virtual displacement. The position can be expressed 
as 

f P = r T e° (2) 

r are the scalar components and e° = [ijk] T is a set of 
unit vectors. Let u be the states of blade deformation. 
The virtual displacement then becomes 

Sfp = (D5u) t e° (3) 

where D is the derivative matrix of the scalar components 
r with respect to the states u. Henceforth, matrices are 
denoted in bold capital letters. In beam theory, the states 
are typically three deflections and three rotations u = 
[iti, i(2, u 3 , 9\, 6*2, 0 3 ]. D is then the (3 x 6) derivative 
matrix. Expressing the force in the same basis 

F = F T e° (4) 

gives 

SW = F T DSu (5) 

If q are the N generalized nodal displacements of a finite 
element containing the point P, and H is the (6 x N ) 
elemental shape function matrix, it follows from Su = 
USq 

SW = F T DUSq = Q T Sq (6) 

The generalized nodal force is then given by 

Q = H t D t F (7) 

The term D 1 transmits the airloads in 3-D space to the 
1-D beam structure. The matrix D varies with the choice 
of beam theory. Note that, from eqn. 3, the virtual dis- 
placement components are 

Sr P = D6u = DH<5 9 (8) 
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Equations 7 and 8 highlight the well-known relation 
that the generalized structural forcing vector relate to 
the aerodynamic forcing via the transpose of the relation 
that connects the aerodynamic deflections to structural 
deflections. Here, DH can be interpreted as the equiva- 
lent elemental shape functions in 3-D space for the cor- 
responding beam theory. 

Note that the force method ensures that the total 
integrated forces (and moments) remain the same when 
transferred from the fluid to the structural domain. How- 
ever, because the point of application of F within each 
surface patch is arbitrary, the method is energy conser- 
vative only in the limit of mesh refinement. Strict con- 
servation can be enforced using traction coupling. 



Figure 2: Fluid pressure and surface shear over a 
rotor blade differential area dA 

The direct traction coupling is formulated as follows. 
Consider a differential area dA within the surface patch 
AH of magnitude cL4 and unit normal n, Fig. 2. 

dA = ndA = n T e°dA (9) 

The differential force generated by the pressure is 

— pdA = — pn T e°dA 

The differential force generated by the fluid stress tensor 
along the direction of dA is 

ctf • dA = of • ndA = (cr-pn) T e° dA 

In the force coupling case, the following integration was 
performed in the fluid domain. 


F= —p dA + up • dA 

Jaa 


/ pn T dA + (c T F n) T dA 

Jaa Jaa 


e° = F T e° 


(10) 


The integrated force was then transferred to the struc- 
tural domain. H and D, required for this calculation, 
were available in the structural domain. Now consider 


the virtual work calculation using fluid traction. 


SW = 


IJ 


—p dA + a-p • dA ) • Srp 


[ [— pn T + (ctfu) t ] e° • (D Su) T e°dA (11) 

Jaa 

f [—pn T + (t7Fn) T ]DH5qdA = Q T 5q 
Jaa 


The generalized nodal force then becomes 


Q= f [-H T D T np + H T D T o- F n] dA (12) 

Jaa 

Note that the integration involves variables from both 
domains. If performed in the fluid domain, exact values 
of D and H must be received from the structural do- 
main at the fluid Gauss points. Similarly if performed in 
the structural domain, exact values of p and a f must be 
received from the fluid domain at the structural Gauss 
points. The value of n is different in both (except in the 
ideal case when the meshes and spatial orders match ex- 
actly). Unlike integrated force coupling, direct traction 
coupling ensures the exact calculation of virtual work. 

There has been significant contributions by vari- 
ous researchers to address the issue of conservation and 
preservation, within the context of temporal accuracy. 
See for example, Maman and Farhat [19], Cebral and 
Lohner [20], Farhat et al. [21], Slone et al [22], and Mich- 
ler et al. [23]. The key conclusion is that exact conser- 
vation and preservation can be ensured in the limit of 
mesh refinement only when interpolations of all the vari- 
ables are performed using schemes consistent with their 
domain. For example, p and ap can be interpolated only 
in the fluid domain and then transferred. Similarly, de- 
formations can be interpolated only in the structural do- 
main and then transferred. 

Neither of these methods have been applied to rotor- 
craft CFD/CSD so far. In the current applications, the 
rotor blade is excited by force and moment distributions 
(per span) along its local sectional elastic axis. The force 
and moment excitations are obtained by integrating the 
fluid stresses along chord- wise strips. Chord- wise strips 
are an easy and natural choice for structured grids. The 
excitations, if calculated at the structural Gauss points 
along the elastic axis (using interpolation in the fluid do- 
main), or imposed as concentrated forces, one for each 
strip using the shape functions, would satisfy conserva- 
tion and preservation requirements in the limit of mesh 
refinement. 

The exact generic methods provide significant ad- 
vantages over the current methods of sectional airloads. 
First, the method of sectional airloads is arbitrary for 
advanced geometry blades, a simple example of which is 
the BERP tip. For such geometries it is necessary to 
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Figure 3: FEM model (MSC NASTRAN) of a SH- 
60 Sea Hawk fuselage developed by Sikorsky and 
used in UMARC rotor-fuselage coupled dynamic 
analysis; Yang et al. 2004 



Figure 4: FEM model of a F-16 fighter used in a 
doublet-lattice aerodynamic model based typical 
fixed-wing flutter analysis; Denegri 2005 


use the exact generic methods for high-fidelity coupling. 
Second, the generic methods are equally applicable for 
both structured and unstructured grids, the latter being 
increasingly applied today for near-blade flows. Third, 
the methods are well-suited for the long term CFD goal 
of near blade adaptive refinement. Fourth, the generic 
methods are equally applicable for the rotor, fuselage, 
and any structure in general, as long as the interface ge- 
ometry representation is clearly defined. 


Interface geometry representation 

The issue of interface geometry representation in- 
volves the rigorous description of the surface based on the 
underlying structural model. The method of description 
differs based on: (1) whether the structural shape func- 



Figure 5: Detailed FEM model of a F-16 fighter 
with 0.17 M DOFs for RANS CFD/CSD based 
aeroelastic response prediction; Farhat et al. 
2003 


tions are available, and (2) even if the shape functions are 
available whether the structural model reaches out to the 
wetted surface. The problem posed by the unavailabil- 
ity of shape functions is a practical complication that is 
ideally avoided. The problem posed by the structure not 
reaching out to the wetted surface is more fundamen- 
tal. For a single component structure, e.g. a rotor blade 
beam model, the necessary interpolation and extrapola- 
tions are rigorously defined by the underlying theory (e.g. 
note the term D in eqn 7 in addition to H). Complica- 
tions arise for multiple intersecting sub-structures, as is 
the case for airframes. 

For illustration, figures 3 and 4 show state-of-the- 
art dynamic models of rotary and fixed wing airframes. 
The details of the work are given later, in the respec- 
tive survey sections. The emphasis in dynamic models 
is on the key load bearing components and internal load 
paths. The emphasis in fluid-structure coupling is on the 
external shell. For such purposes models which include 
the outer shell, and in addition, the internal load paths, 
are most appropriate. For example, Fig. 5, shows a more 
detailed model of the same F-16, including the outer shell 
(details of the work is cited later). 

Detailed modeling of the fuselage for fluid-structure 
interaction has not been a key focus so far in rotorcraft 
for the following reason. The dominant dynamics of a 
helicopter fuselage, unlike fixed-wing, is vibration, which 
occurs at high frequencies (usually pNf,/ rev). The domi- 
nant source of this vibration is not the fuselage flow field, 
but the shaft transmitted forcing from the rotor. The 
forcing from the rotor primarily stems from the rotor 
flow field, the interactional effect of the fuselage on the 
rotor flow field occurs, nominally, at lower frequencies 
(1/rev). 

There are important conditions where the above ar- 
gument is invalid. First, at low speed flight (around 
p = 0.15) direct impingement of the rotor tip vortices 
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(a) 1-D line fuselage 



(b) 3-D detailed fuselage 

Figure 6: Measured and NASTRAN calculated 
fuselage natural frequencies of the AH-1G heli- 
copter; Yeo and Chopra 2002 


Figure 7: FEM line and 3-D fuselage NAS- 

TRAN models of the AH-1G helicopter from the 
DAMVIBS program used for rotor-fuselage cou- 
pled dynamic analysis; Yeo and Chopra 2002 


on the tail boom is a key mechanism of vibration (at 
Nb/ rev). RANS was unable to resolve this flow field thus 
far, but is beginning to do so now. Second, for advanced 
configurations with low shaft clearance (for low drag) di- 
rect wake impingement is a significant source of vibration 
for a wide range of flight speeds. There has been limited 
resolution of the interactional flow field around the fuse- 
lage thus far. RANS based CFD, and alternative wake 
CFD models like the Vorticity Transport Model [24, 25], 
are beginning to do so. The emergence of interactional 
CFD opens opportunity to calculate vibration and buffet 
loads in response to the interactional flow field. 

Reproducing the structural frequencies of the fuse- 
lage (need at least up to 30 Hz) is in itself difficult in 
helicopters (more than 20% error without detail model- 
ing) . Difficult components lie in critical load paths (main 
rotor hub/pylon/transmission case). Dynamicists often 


use line fuselages to re-produce measured frequencies. 
For example, the calculated frequencies of an AH-1G 
(two-bladed, teetering) fuselage for two different types 
of models are compared with measurements in Fig. 6. 
The corresponding models are shown in Fig. 7. The sim- 
pler model is designed to generate similar frequencies as 
the more detailed model, and is often preferred. The 
fluid-structure coupling methodology, however, is more 
involved for the simpler model compared to the detailed 
model. The need for extrapolation, a key source of error, 
is minimized in the case of detailed models that include 
the external shell structure. Models with both internal 
and external details of the structure are used regularly by 
crash researchers, for very short-time high impact tran- 
sients (less than 0.1 sec). An example of a crash test 
simulation is shown in Fig. 8, taken from Ref. [26]. 

Constructing a generalized modular interface which 




(a) Detailed shell fuselage 


Acceleration, g 



(b) Measured and predicted acceler- 
ations at a front left fuselage station 

Figure 8: Detailed FEM composite fuselage (MSC 
Dytran) as used in high-fidelity crash models; 
Fasanella et al. 2007 (courtesy Karen Jackson, 
NASA Langley) 


is adaptable to alternative levels of underlying struc- 
tural fidelity requires innovative interpolation and ex- 
trapolation methods for interface displacement mapping. 
A large volume of literature has been devoted to such 
methods for fluid-structure coupling. Broadly, they are 
divided into two categories. First, the surface tracking 
methods. These rely on the underlying shape functions 
while trying to mitigate the complexities associated with 
multi-component sub-structures. Second, the surface fit- 
ting methods. These are focused more on independent 
(from fluid and structural solvers) formulations without 
the need for shape functions. In fixed wing applications, 
the current practice is to use a combination of these in- 
terpolation methods to map deformations from CSD to 


CFD. Mapping airloads from CFD to CSD then simply 
amounts to the correct evaluation of virtual work based 
on the deformation mapping (compare eqn. 8 with eqns. 7 
and 12). 

A displacement mapping must: (1) recover large de- 
flections and rotations of the underlying structure from 
the surface displacement, and (2) ensure displacement as- 
sociation with time, i.e. attachment points must remain 
attached throughout the simulation. Arbitrary proce- 
dures like a nearest neighbor mapping, although simple 
and fast, do not meet these requirements. Moreover, it 
prevents the correct evaluation of the virtual work term. 
Procedures which rely on underlying structural shape 
functions for multi-component blending and extrapola- 
tion meet these requirements. They are called surface 
tracking methods as mentioned earlier. When the shape 
functions are not available, the most general procedure 
is to use radial basis functions. These are also called 
surface fitting methods. The well-known surface fit- 
ting methods include the Multiquadric-Biharmonic (MQ) 
method, the method of surface splines commonly called 
the Infinite Plate Spline (IPS), the method of Thin 
Plate Spline (TPS), and Non-Uniform Rational B-Spline 
(NURBS). These methods have been systematically stud- 
ied by Hounjet and Meijer [27] in 1994, and Smith et 
al. [28] in 2000 (see references therein for details of the 
above methods). The study by Smith et al. found TPS 
to be most suitable for interpolation and extrapolation 
of generalized structural geometries. These radial ba- 
sis functions are defined on the whole domain of a sub- 
structure. Complications arise for multiple intersecting 
sub-structures in 3-D space. 

Fitting methods applicable in 3-D, which do not use 
radial basis functions, have also been implemented, no- 
table among which are the Constant Volume Tetrahe- 
dral (CVT) method by Goura et al. [29, 30], and the 
Inverse Boundary Element (IBE) method by Chen and 
Jadic [31]. Even though applicable in 3-D, these methods 
are less suitable for large deformation nonlinear prob- 
lems. The first method provides a nonlinear mapping. 
Linearization destroys the exactness of rigid body rota- 
tions. The second method requires iterations to conform 
to the deformed geometry, unless the deformations are 
assumed to lie within the range of linear elasticity. 

An emerging technique is the use of 3-D radial basis 
functions with compact support. These are defined lo- 
cally, and hence are easily applicable for surfaces with 
high 3-D curvature. The method is based on multi- 
variate scattered data interpolation. Originally proposed 
by Wu [32] and Wendland [33] , and subsequently applied 
by Beckert and Wendland [34] for fluid-structure inter- 
action problems, the method has found increasing appli- 
cation in recent large-scale fixed-wing CFD/CSD calcu- 
lations [35, 36]. 
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FIXED WING AIRCRAFT 

A comprehensive survey of the state-of-the-art in 
fixed-wing aeroelasticity can be found in Livne [37] in 
2003. Schuster et al. [38], in the same year, have reviewed 
the state-of-the-art in computational aeroelasticity. The 
term was used broadly, covering unsteady aerodynam- 
ics to CFD and beam models to detailed FEM. In 2004, 
Kamakoti and Shyy [39] described the RANS/FEM cou- 
pling procedures used for fixed wing applications. The 
AGARD 445.6 wing as taken as the baseline configura- 
tion for illustration. More recent developments in com- 
putational aeroelasticity are highlighted by Bartels and 
Sayma [40] in 2007. In the same year, a report detailing 
European research towards the development, evaluation, 
and transition of high-fidelity aeroelastic tools for full 
aircraft has also been published [41]. 

The following sections review the status of compu- 
tational aeroelasticity in fixed wing applications. The 
first section briefly describes the status of detailed CSD. 
The second section summarizes the important aeroelas- 
tic phenomena and the status of high-fidelity CFD/CSD 
methods in predicting them. 

Structural Modeling 

First, the meaning of ‘detailed structures’ must be 
clarified. Shown in Fig. 9 is a typical Boeing 737 wing sec- 
tion with a simple main/aft double flap system with a sin- 
gle slotted thrust gate, in its retracted position [42] . Only 
the trailing-edge is shown, the three position Kruger slats 
in the leading edge are not shown. Typical structures 
such as this, incorporate kinematic constraints to rigid 
body modes to flexible deformations. For aeroelastic and 
aeroservoelastic studies, a distinction is made between 
global and local behavior of structures. Only the global 
model is coupled to the external flow. The loads gen- 
erated by the global model are then transferred to the 
internal structures wherever needed, locally. 

The local models handle details, which include FEM 
of individual composite ply-layups, rigid body mecha- 
nisms, redundant load paths, local buckling, crack prop- 
agation, stress concentration due to manufacturing im- 
perfections, embedded cooling, and integrated smart ma- 
terials and actuators. Local effects (e.g. actuator loads, 
local buckling, composite tailoring) are included via lower 
order models. The advent of HPC opens opportunity to 
couple the important frequencies of local analysis directly 
to global analysis using detailed modeling. Performed ju- 
diciously, based on a fundamental understanding of the 
key phenomena, it can lead to enormous payoffs in design 
innovation and cost. 

A comprehensive treatise on fixed-wing structural 
analysis and design, as performed in the Boeing Corn- 



Figure 9: The outboard wing cross section of a 
Boeing 737-NG aircraft; van Dam 2002 



Figure 10: A Boeing 777 FEM stress model with 
details of major internal components; Mohaghegh 
2005 

pany can be found in Mohaghegh [43]. The philosophy 
of fail safe design has led to the well-known practice of 
discretely stiffened wing and fuselage panels and mul- 
tiply redundant two and three-piece primary bulkheads 
in the construction of modern aircraft structures. To- 
day, global analysis can be performed using full aircraft 
FEM and RANS. The structural model involves large- 
scale but usually linear FEM models. The global loads 
are then used for internal stress calculations using de- 
tailed FEM of the type shown in Fig. 10. Most often, 
this step involves a static analysis. Since the 1990s, all 
primary structures of commercial airplanes like the B777 
and A340 are certified using such FEM analysis. A typ- 
ical example is the wing-body junction which includes 
the cargo bay and landing gear — a region of complex 
redundant load paths and structural discontinuities, see 
Fig. 11. The need for analysing such structures led to the 
original development of sub-structuring methods. Com- 
posites materials are used on selected components like 
the radomes, fairings, all control surface panels, engine 
nacelles, and torque boxes. The main purpose is to pro- 
vide higher impact, fatigue, and corrosion resistance. 

Almost every modern rotor blade is an all compos- 
ite construction (although the main spar is frequently 
of Titanium construction). The internal layout of rotor 
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Figure 11: A routine Boeing 747 wing-body interaction FEM stress analysis model with typical sub- 
structures; Mohaghegh 2005 


blades, however, are much simpler compared to the me- 
chanical complexity of fixed wing sections. The control 
mechanism for conventional main rotors are all essen- 
tially concentrated at the root end. In the case of rotors, 
the state-of-the-art in global and local models are ID 
beam FEM and 3D stress analysis respectively. The ID 
beam properties are first obtained by a separation of the 
3D problem into a ID problem and a 2D cross-sectional 
analysis problem (see section on rotary wing) . Analogous 
to its aerodynamic lifting-line counterpart, the theory is 
not meant for end stresses or stress concentrations at 
the attachment point. More importantly, the root end 
flex-beams and torque tube structures of modern hin- 
geless and bearingless rotors require that the external 
loads be transferred correctly to internal stresses. This 
is currently performed by a detailed re-calculation using 
commercial FEM codes like ABAQUS, ANSYS, etc, in a 
static fashion (using the dynamic forcing obtained from 
the global model) . It is desirable that future HPC based 
simulations make this approximate re-calculation redun- 
dant. 

While the 1-D beam FEM model of a rotor blade 
is comparatively simpler, the models include geometric 
nonlinearities rigorously that are critical for rotors. Ro- 
tary wing aeroelasticity is inherently nonlinear, the non- 
linearities stem from the centrifugal and Coriolis forces 
due to blade rotation. The Coriolis nonlinearity is in turn 
coupled with the rotor trim state. This level of coupling 
between CSD and vehicle dynamics is absent in fixed 
wings, and this has its implications on the efficiency of 
methods that couple rotorcraft CSD with CFD. 

Aeroelastic Phenomena and CFD/CSD Predic- 
tions 

Flutter 

Flutter analysis in subsonic regimes are satisfacto- 
rily performed by standard k and p — k methods using 
doublet-lattice type linear aerodynamic models. These 
options are available today as part of FEM codes like 
NASTRAN and ASTROS. In supersonic regimes, meth- 


ods related to piston theory are used. A description of the 
current status of these methods, along with their aeroe- 
lastic applications, can be found in Yurkovich [44]. High- 
performance fighter aircrafts are, however, flutter critical 
in the transonic regime where unsteady shock motions 
of three different types determine the nature of energy 
transfer between the fluid and the structure. A linear 
structure is still enough, only the aerodynamic nonlin- 
earities need to be predicted. Hence high-fidelity RANS 
CFD approaches are sought. 

Unlike helicopter blades, where flap-lag flutter of 
hingeless rotors is determined by structural nonlineari- 
ties coupled to the vehicle operating state, fixed wing 
bending-torsion flutter is a linear phenomena. The effect 
of vehicle g-level is only via the effect of trim angle of at- 
tack on the shock motion. At a given angle of attack, and 
dynamic pressure, a conceptually straight-forward fluid- 
structure transient response can be used to extract the 
frequency and damping. Because classical logarithmic 
decay methods are sensitive to noise, refined eigensys- 
tem realization algorithms are used to extract these pa- 
rameters for the particular modes of interest. Most often 
(except for free-body control surface flutter) these are the 
first few antisymmetric modes of the aircraft. Depend- 
ing on the convergence and divergence of response, gen- 
erally 5 to 6 transient calculations are enough to bracket 
a flutter point (i.e. a frequency, damping vs. dynamic 
pressure point corresponding to zero damping at a given 
Mach number). 

One of the early CFD based flutter calculation for 
a complete aircraft was performed by Melville in 2001 
and 2002 using Euler [45] and RANS [46] coupling, for a 
F-16 fighter aircraft. The linear FEM structural model 
used was the same as used by Denegri earlier [47] for a 
doublet-lattice based flutter analysis. Ten symmetric and 
ten anti-symmetric mode shapes were considered. An ex- 
tra mode was incorporated for the leading edge flap. The 
flap was represented as a damped first-order system that 
relaxed into a commanded position. The FEM model 
is shown in Fig. 12(a). It was built from several sim- 
pler models representing each of the main components 
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(a) FEM structural model 



(b) Fluid-structure interface ge- 
ometry representation 



(c) Measured and predicted flutter boundary 


Figure 12: Measured and predicted flutter onset boundaries for a F-16 fighter aircraft using 

RANS/FEM modal coupling; Angle of attack a=1.5°; Melville 2001, 2002 



(a) Predicted and measured torsion mode damping 



(b) Predicted and measured bending / torsion frequencies 


Figure 13: Predicted and measured aeroelastic parameters for a F-16 fighter aircraft at 3000 m altitude 
using RANS/full FEM coupling; Farhat et al. 2003 


the fuselage, wings, leading-edge flap, flaperon, tip 
launcher, and tails. The fluid-structure interface repre- 
sentation was constructed by seperate extrapolation of 
each component. The interface is shown in Fig. 12(b). 
The curvature was provided by the underlying structural 
components, extrapolation to surface was done linearly. 
The RANS domain (Baldwin-Lomax) contained 20 to 36 
overset grids with 1.9M to 3.5M points. The details of 
the inlet, exhaust, ventral fins and under-wing store were 
not included. The coupling method and the moving grid 
mechanism (with its associated Geometric Conservation 


Law) was based on an earlier work by Morton et al. [48] . 
A single CPU was used for structures and grid deforma- 
tion, the fluids solver was run in parallel. 

A perturbation response damped the three primary 
symmetric modes into steady deflections. The three 
primary anti-symmetric modes were lowly damped and 
showed zero mean oscillations. The conventional proce- 
dure for extracting modal damping is to use classical log- 
arithmic decay methods. For a full FEM model, real-time 
parameter identification methods, as used in flight tests, 
are used. The procedure is recognized as a challenge in 
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(a) A Boeing twin-engine transport wind-tunnel model 
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(b) Transonic flutter onset dynamic pressures for a hump mode and a tip mode 


Figure 14: Predicted and measured flutter onset for a twin-engine Boeing transport aircraft model 
using RANS/FEM modal coupling; Hong et al. 2003, reproduced from Bartels and Sayma 2007 


general. Often, a doublet-lattice based linear analysis 
is performed a priori to identify the critical modes. In 
this particular instance for example, linear analysis re- 
sults from Denegri [47] was used to identify the second 
mode as the critical mode. This mode was then tracked 
during the CFD/CSD simulation. Figure 12(c) compares 
the predicted flutter boundary with flight test data. 

Farhat and his co-researchers [49, 50] have per- 
formed RANS/full FEM simulations of the same aircraft. 
The structural model was that shown earlier in Fig. 5. It 
contained around 0.17M linear DOFs comprised of bar, 
beam, solid, plate, shell, and composite elements. The 
fluid solver was run on up to 24 CPUs on an SGI Ori- 
gin 3200 computer with parallel speed up of 21. Of the 
total CPU time, 60% was spent on fluids, 38% in mesh 
deformation, and only 2% on structures. The authors 
note that for nonlinear (geometric) FEM, the structures 
require 15% of time. The total time for one transient 
calculation on 24 CPUs was 2.5 hours. An Eigensystem 
Realization algorithm was used to extract the frequencies 


and damping of the two lowest modes, see Fig. 13. The 
algorithm requires only two cycles of history as long as 
the sampling rate is 500-1000 Hz. A typical wing-section 
analysis using 2-D CFD coupled to a 2-DOF spring-mass 
system was also performed for comparison. Recently in 
2007, Lieu and Farhat [51] have used the same test case to 
systematically compare these results with those obtained 
using Proper Orthogonal Decomposition based reduced 
order aerodynamic models. 

A description of the status of RANS based flutter 
predictions in commercial aircraft design can be found in 
the review paper by Bartels and Sayma [40] . A collabora- 
tive Boeing-NASA effort recently examined flutter onset 
on a half-span twin-engine transport model including the 
aerodynamic interactions of the wing, nacelle, and strut, 
the effect of wind tunnel wall, and the effect of static aero. 
The details of the work originally appeared in Hong et 
al [52], and are reproduced here in Fig. 14 from refer- 
ence [40]. The CFD code used was CFL3D (version 6), 
which coupled RANS to structural modes. A second- 
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(a) Flight test 


(b) 1/72-scale model test 


Figure 15: Flow visualization of LEX vortices; (a) An LEX vortex on a F/A-18 High Angle-of-attack 
Research Vehicle (HARV), wing-tip view at a=25°, /?=-1.4°; (b) A 1/72-scale model plan view at 
o=30°; Lee 2000 


order predictor-corrector based coupling was used in this 
study between the fluid and structure. The model exhib- 
ited a ‘hump’ mode, i.e. a lightly damped mode which 
is unstable only on a limited range of dynamic pressure. 
The analysis predicted this mode. Predictions for an- 
other unstable mode, the tip mode, showed inconsistent 
results. Predictions showed greater error when the ac- 
tual cruise shape was used compared to the undeformed 
jig shape. 

Tail Buffet 

Fixed-wing tail buffet studies date back to the 1930s 
and officially began with a Junkers F13 ge monoplane 
crash in Meopham in England. The official British inves- 
tigation attributed the tragic accident to horizontal tail 
buffeting. Subsequently, early work focused mostly on 
horizontal tail buffet. Since the late 1970s and 80s the 
buffet problem has received renewed attention with the 
introduction of the F-14, F-15, F-16 and F/A-18 fighters 
in the US Air Force and Navy. This time, the focus has 
been on vertical tail buffet. 

Modern fighters are designed with maneuver capa- 
bilities at high angle of attack (AOA), typically 25° to 
35°, sometimes as high as 60°. To generate lift at these 
high AOA, a Leading edge Extension (LEX) is designed 
on each wing to create two vortices that trail downstream 
over the aircraft. At the same time small disturbances 
from the fore-body are also reinforced by the two LEX 
vortices. The vertical tails (also called fins) which pro- 
vide directional stability and control, are placed to take 
full advantage of these two vortices. However, if the vor- 


tices fail to traverse the adverse pressure gradient on the 
wing, and burst, they impinge on the tail instead and 
produce large vibration and consequent fatigue damage. 
This phenomena is called tail buffet. Figure 15 shows 
flow visualization of vortical flows that produce tail buf- 
fet. 

Tail buffet without the interaction with LEX vortex 
have also been reported. For example, on the F-15, the 
main wing vortex system generated buffet loads. It is 
only for the F/A-18 aircraft that a large volume of ex- 
perimental and analytical studies exist for a systematic 
understanding of the problem. Lee [54] in 2000 has com- 
prehensively described the problem, and reviewed the 
status of fundamental understanding, experimental data, 
and analytical predictions of the tail buffet. The problem 
has acquired increased relevance as the three largest U.S. 
fighter aircraft programs have twin tails, F/A-18E/F, F- 
22, and the F-35 (Joint Strike Fighter). 

There have been several noteworthy attempts at al- 
leviating tail buffet from the late 80s through 90s using 
extensive wind tunnel and flight tests. The only method 
used today is the LEX fence developed by McDonnell 
Douglas Corporation via trial and error using wind tun- 
nel experiments and flight tests. The LEX fence is a 
flat trapezoidal plate standing on the LEX upper surface 
aligned in a stream-wise direction. 

The complexity of the problem and its sensitivity 
to flight conditions and aircraft configurations makes it 
ideally suited for a high-fidelity analysis. At the time of 
Lee’s review there was only limited RANS calculation of 
the tail buffet problem. These early calculations were by 
Rizk and Gee [55] in 1992, Ghaffari et al [56] in 1993, and 
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(a) Instantaneous streamlines of LEX and wing 
vortices without LEX fence 


(b) Instantaneous streamlines of LEX and wing 
vortices with LEX fence 
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Figure 16: Predicted and measured tail bending moments during F/A-18 tail buffet alleviation us- 
ing LEX fences; RANS/FEM coupling, FEM of vertical tail only; Bending moments are integrated 
differential pressure moments, not structural gauge loads; Sheta 2004 


Gee et al. [57] in 1996. These studies focussed on CFD, 
the structure was assumed to be rigid, and had limited 
success in predicting structural buffet loads on the tail. 

A RANS/FEM prediction of buffet loads on the 
F/A-18 tail was carried out by Sheta in 2004 [58]. The 
study clarified the fundamental contribution of the LEX 
fence in alleviating buffet loads. The aircraft was as- 
sumed to be rigid, but the tail was modeled using 3- 
D hexahedral brick FEM (isotropic, Aluminum). The 
RANS domain had 53 structured blocks with a total of 
2.68 M grid points. The total time required on 6 CPUs 
(of 1.0 GHz speed each) was 142 h. Figure 16 shows the 
predicted flowfields and tail bending moments from this 
study with and without the LEX fences. The off-body 
flow was considered to be laminar in this study. 

Recently, Morton et al. [59] have used detached- 


eddy simulation (DES) based turbulence model for a 
more meaningful calculation of vortex breakdown. DES 
is a popular RANS/LES hybrid model, ft uses RANS 
in the boundary layer and LES in the regions of mas- 
sive separation. Designed primarily for external flows, 
it can be implemented by a relatively simple modifica- 
tion of an existing Spalart-Allmaras (S-A) RANS model 
in an unsteady solver. The LES domain does not re- 
quire an explicitly defined SGS model. The modification 
to the S-A model automatically implies the existence of 
a Smagorinsky-type SGS model under the assumption of 
local equilibrium between the production and dissipation 
of eddy viscosity. The simulations were performed for the 
F-18 High Angle-of- Attack Research Vehicle (HARV). 
For the chosen flight condition, the F-18 HARV configu- 
ration matches an F/A-18C configuration with leading- 
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(a) LEX vortex breakdown around a 
F/A-18 aircraft using DES model 
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(b) Vortex breakdown location (nose is x=0) 
validated with flight test data 


Figure 17: RANS calculation of vortex bursting 
and tail impingement using Detached Eddy Sim- 
ulation (DES); no LEX fences; Morton et al. 2007 


edge flaps set to -33°, trailing-edge flaps set to 0°, and 
the diverter slot through the LEX fence closed. The an- 
gle of attack was 30°, Mach number 0.2755, Reynolds 
number 13M (corresponds to altitude of 20,000 ft). The 
DES simulations were performed on a baseline unstruc- 
tured grid of 3.6M cells (using the commercial CFD code 
Cobalt). A time-averaged solution was used to produce 
an AMR grid with 3.9M cells using Pirzadeh’s grid adap- 
tation method [60]. The analysis, however, did not in- 
clude a structural dynamic model; the focus was on the 
flow field. 

Limit Cycle Oscillations (LCO) 

Fundamental understanding of LCO and its satis- 
factory prediction remain beyond the state of the art in 
fixed-wing aeroelasticity. The LCO is a general term. In 


fixed wing fighters it is characterized by sustained, pe- 
riodic, antisymmetric oscillations of the wing which are 
not catastrophically divergent. Inside the fuselage, it is 
felt as a lateral motion. It has been a persistent prob- 
lem since the mid-1970s in high subsonic and transonic 
speeds for configurations with missiles on the wingtips 
and heavy stores on the outboard pylons. The F-16 and 
F/A-18 fighters with wing tip missile launchers tradition- 
ally encountered LCO. Similar phenomena are expected 
for the F-22 and F-35 as flight tests begin. 

Bunton and Denegri [61] presents an insightful com- 
mentary on fixed-wing LCO. The phenomenon is similar 
to classical flutter, in that the oscillations occur at a sin- 
gle frequency close to that predicted by linear flutter. 
The mode shape, at the onset, also resembles that pre- 
dicted by linear flutter. In level flight, the onset speed 
lie quite close to that predicted by linear flutter. As a re- 
sult LCO has been termed limit cycle flutter or constant 
amplitude flutter by many, a terminology not preferred 
by others as the word ‘flutter’ has historically referred to 
catastrophic divergence. Many speculate that the gene- 
sis of LCO is not the classical bending-torsion flutter at 
all, but the aircraft flight control system. 

During LCO, the amplitude remains constant for a 
given, stable flight condition. If the aircraft is accelerated 
to a higher speed, the amplitude would increase and set- 
tle into a higher level once the higher speed is reached. If 
an attempt is made to back out of LCO, the oscillations 
often persist below the original onset values (hysteresis) . 
In addition to LCOs involving the entire aircraft-wing 
store structure, LCOs can occur involving localized non- 
linearities like control surface free-play. 

LCO occurs both in level flight as well as in high 
load factor maneuvers. Usually the onset speed is lower 
in the latter. Compared to LCO onset in level flight, 
those occurring in maneuvers are much less predictable. 
For example, LCO magnitudes can increase from lg to 4g 
and then begin decreasing so as to disappear at 5.5g [61]. 
For a different store configuration there may not be any 
LCO at all until 5g, at which LCO begins, and grow 
so rapidly so as to resemble the onset of classical catas- 
trophic flutter. The key unresolved challenge, beyond 
that of predicting the onset, is predicting the magnitude 
of LCO. In fixed wing, unlike rotorcraft, the fatigue of 
components is of lesser concern compared to the limit 
loads. (Fatigue concerns are evident however for test air- 
craft repeatedly exposed to flutter tests and LCO.) LCO 
concerns are related to weapons accuracy for smart and 
unguided weapons, and reliability of ordnance release. 
The transonic aerodynamic nonlinearity that can lead to 
LCO has been explained by Bendiksen [62] based on the 
categorization of shock motions by Tij deman and Seebass 
in 1980 [63]. It was shown that, for a clean wing, when 
a bending-torsion flutter mode triggers a transition from 
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(a) External stores on a F-16 fighter during flutter and LCO 
flight tests 
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(b) Measured maximum wing-tip LCO response in 
level flight 

Figure 18: Measured wing-tip LCO deformations 
of a F-16 fighter aircraft with under-wing stores; 
Denegri et al. 2005 


(a) Configuration 1 



(b) Configuration 2 

Figure 19: Predicted and measured LCO magni- 
tudes using RANS CFD/FEM coupling in level 
flight; (a) Effect of slendor body-wing model for 
launchers, (b) LCO predicted only with launcher 
model; Thomas et al. 2007 


type-A (shock motion oscillates sinusoidally) to type-B 
(shock disappears during a downstream motion) shock 
motion, then the associated rapid decrease of aerody- 
namic work done per cycle can lead to LCO. The role 
played by under-wing stores and structural nonlineari- 
ties including free-play are more configuration specific. 
As flight testing becomes more and more expensive, and 
the number and combinations of under-wing store con- 
figurations increase, higher-fidelity methods are ideally 
suited for a comprehensive attack on the LCO problem. 

One of the early high-fidelity attempt to simulate 
LCO was that of Melville in 2002 [46] described earlier 
in the section on flutter. While calculating the onset of 
flutter he investigated whether the undamped response 
(or lowly damped response near the flutter boundary) 
resembled an LCO as seen in flight. What distinguishes 
an LCO response from a zero-damping periodic flutter 


point is that the LCO amplitude is independent of the 
initial excitation, and depends only on the underlying 
limiting mechanism. The second mode excitation, when 
calculated with different set of initial excitation, settled 
to a different amplitude (proportional to the initial exci- 
tation). Thus it was not LCO. 

More recently, in 2005, Denegri et al. [64] have gener- 
ated benchmark LCO data for the F-16 aircraft. A sam- 
ple of benchmark measured LCO deflections along with 
the associated wing measurement locations are shown in 
Fig. 18. Oscillation frequencies range from 7.9Hz to 8.2 
Hz with lower frequencies around Mach 0.9. 

In 2007, Thomas et al. [65] have used a nonlinear 
harmonic balance based RANS solver to calculate the 
LCO magnitudes for different configurations of the F-16. 
A slender-body wing-theory was used to model the ex- 
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ternal weapons and wing stores and study their impact 
on LCO prediction for eight different weapons and store 
configurations. A comparison of the predicted LCO mag- 
nitudes with those of measured data, as in Fig. 18(b), is 
given in Fig. 19(a). It corresponds to 2000 ft altitude, 
1.5° angle of attack flight test. The figure is sufficient 
to highlight the poor state of affairs in LCO prediction. 
Figure 19(b) compares predictions for a different config- 
uration, but at the same altitude and angle of attack. 
The LCO could at least be predicted with the wing store 
model. 

Parker et al. [66] have carried out systematic studies 
on the effect of under-wing stores on LCO using a simple 
rectangular cantilevered wing model. The details of the 
wing model can be found in [67]. The model is based 
on Eastep and Olsen’s heavy version [68] of the origi- 
nal Goland wing [69]. They found that the role of store 
aerodynamics was to transfer additional energy into the 
structure increasing the magnitude of limit cycle oscilla- 
tion. 

In addition to the F-16 and F/A-18, LCO studies 
have been documented on a B-l like configuration [71] 
and on a B-2 bomber. The latter phenomena is now well- 
known as the Residual Pitch Oscillation (RPO) and has 
been studied by Dreim and his co-researchers [72] using 
the Computational Aeroelasticity Program - Transonic 
Small Disturbance (CAP-TSDv) code. 



Figure 20: FEM model of bladed disks: (a) a disk 
sector; (b) blade contact surface; (c) disk contact 
surface; and (d) an area contact element; Petrov 
and Ewins 2006. 


TURBOMACHINES 

The advances in fundamental understanding, analy- 
sis, and prediction of turbomachinery aeroelasticity can 
be broadly traced from 1970s to the present time by sur- 
vey papers in Mikolajczak et al. in 1975 [73], Bendikson 
in 1990 [74], Marshall and Imregun in 1996 [75], and Bar- 
tels and Sayma in 2007 [40]. The historical evolution of 
the aerodynamic methods has been traced recently by 
Cumpsty and Greitzer [76]. Structural dynamics (and 
flutter) has been surveyed by Srinivasan in 1997 [77], 
Slater in 1999 [78], and more recently Castanier and 
Pierre in 2006 [79]. 

The purpose of the following sections is to review 
the status of turbomachinery high-fidelity computational 
aeroelasticity. As in the case of fixed wing, the focus here 
is again on CSD, and coupled high-fidelity CFD/CSD 
methods. The first section summarizes the status of 
detailed CSD. The second section summarizes the key 
aeroelastic phenomena, and the status of CFD/CSD in 
predicting them. 

Structural Modeling 

Finite element models of turbine blades comprise of 
millions of DOFs made up to 3-D linear solid elements. 
Commercial FEM packages like ANSYS and ABAQUS 
are routinely used for industrial as well as research pur- 
poses. Turbine blades are made out of nickel based single- 
crystal super alloys. In addition to unsteady stresses, 
turbine blades have to endure high temperature (par- 
ticularly during start-up and shutdown) and exposure 
to high-pressure hydrogen (in case of internal coolant). 
Furthermore, the blades are prone to severe fretting by 
the root end friction damper, and contact damage at the 
dovetail attachments which can lead to crystallographic 
initiation and crack growth. Detailed 3D FEM analyses 
are carried out with special surface-to-surface contact el- 
ements for such purposes, see for example Arakere et 
al. [80] and the references therein. 

The forced response of bladed disks are affected by 
contact surfaces at the blade-disk dovetail attachments, 
part-span interlocking shrouds, and under-platform fric- 
tion dampers at the root. These components are used to 
avoid dangerous resonance regimes. Significant amount 
of research is conducted to model these discontinuity phe- 
nomena within 3-D FEM. Recent work on 3-D FEM mod- 
eling of dovetail attachments can be found in Beisheim 
and Sinclair [81]. Recent work on modeling of under- 
platform friction dampers can be found in Petrov et 
al. [82, 83]. Figure 20, taken from the latter work, illus- 
trates the placement of an area contact element within a 
typical FEM model of a bladed disk. 

When all sectors of a bladed disk are similar, the 
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(a) Three nodal diameter mode of a tuned bladed disk 


(b) Localized mode shape of a mistuned bladed disk 


Figure 21: Mode shapes of a tuned and mistuned modern bladed disk; Castanier and Pierre 2006 


structure is called tuned, or having cyclic symmetry. The 
analysis can then be simplified to one or a few sections. 
Mistiming destroys cyclic symmetry and produces mode 
localization in a bladed disk. This phenomena is partic- 
ularly important for High Cycle Fatigue (HCF), where 
high localized stresses can lead to failure. Figure 21 
shows an example of a modern bladed disk FEM model 
showing the effect of mistiming on mode shapes. Studies 
on mistiming have received attention since 1970s, led by 
the pioneering work of Whitehead [84] and Ewins [85] . 
However, it is only recently that the emergence of whole 
annulus multi-row CFD have made reliable forced re- 
sponse predictions of mistuned bladed disks appear fea- 
sible in near future. For a recent review of the vari- 
ous modeling and analysis techniques, and the assess- 
ment of mistuning sensitivity, see Ref. [79]. The mode 
localization due to mistuning, as mentioned earlier, is 
more pronounced in the presence of frequency veering. 
Frequency veering is a common phenomenon in rotat- 
ing structures. Coupled modes approach each other with 
increase in RPM and veer away after exchanging mode 
shapes due to changing contributions of centrifugal load- 
ing in each. The same phenomena occurs in helicopter 
blades, except that the low aspect ratio of the turbine 
blades (around 1.5-3 compared to 10-15 for rotorcraft) 
and the large flexible disk compress a larger number of 
vibration modes within a smaller frequency range which 
are richer in interaction. A unique phenomena, not to 
be found in rotorcraft blades, is the interaction between 
disk dominated and blade dominated modes. It is ob- 


served, that the sensitivity of a bladed disk to mistuning 
increases with frequency veering. Recent work on mis- 
tuning and mode localization can be found in Petrov and 
Ewins [86], Rivas-Guerra and Mignolet [87], and Kenyon 
et al. [88, 89]. 

The whole annulus treatment of a bladed disk, with 
or without mistuning, is in itself an idealization which 
holds true only for isolated, dynamically independent 
rotors. A disk with identical blades is not necessar- 
ily tuned, particularly for multi-stage compressor blades. 
The modes of one stage interact with those of its neigh- 
bor, which usually contains a different number of blades 
giving rise to mistuning of the coupled stages. Unlike 
helicopters (co-axial) the stage to stage spacing is only a 
fraction of the blade chord. In addition, the disk is large 
and flexible, compared to the hub in helicopters. Bladh 
et al. [90] have studied the effects of structural coupling 
of multi-staged bladed disks, see Fig. 22. 

Aeroelastic Phenomena and CFD/CSD Predic- 
tions 

Aeroelastic Phenomena 

A modern turbofan engine and a schematic illustrat- 
ing its instability prone components are shown in Fig. 23. 

Static aeroelasticity in turbomachines deals with the 
calculation of blade un-twist and un-camber because of 
high centrifugal and pressure loads. There is no diver- 
gence phenomenon. Typical un- twists can be around 5° 
and changes the stagger angle. The manufactured shape 
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Stage 1 Multi-stage Stage 2 

12 blades 16 blades 




(a) FEM for single and two-stage rotor models ( b ) Stage 1 forced response from Engine Order 10 excitation, 

using tuned and mistuned FEM single and two- stage models 


Figure 22: Effect of multistage rotor structural coupling in turbomachines; Bladh et al. 2003 
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(a) The GE 90-115 B growth engine of the Boeing 777 


(b) Schematic of a modern turbofan 


Figure 23: The components of a modern, high bypass ratio turbofan engine 


(cold) is calculated based on the anticipated un-twist and 
un-camber of the operating shape (running). The effect 
is pre-dominant in fan blades. A recent study can be 
found in Wilson et al. [91]. 

Forced response is most critical for turbine blades. 
The High Cycle Fatigue phenomenon is a major issue 
is gas turbine blades. In addition to unsteady forcing, 
high thermal stresses and internal coolant pressure lead 
to fatigue failure. The U.S. Air Force estimates that 55% 
of jet engine safety Class A mishaps (over $1 million in 


damage or loss of aircraft) and 30% of all jet engine main- 
tenance costs are due to HCF [92]. The forced response 
phenomenon also occurs in compressor blades. But here, 
it occurs not only at Engine Order (EO) excitations (sim- 
ilar to blade passage frequency in helicopter main rotors), 
but also at Low Engine Order excitations (LEO). The lat- 
ter arise due to stage-to-stage interactions. Compared to 
turbines, compressors typically have greater number of 
blade stages which produce a greater number of engine 
order excitation frequencies. 
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Figure 24: A typical compressor map showing the 
various flutter boundaries 


Flutter occurs in fan blades, the front stages of axial 
compressor blades, and occasionally in the low pressure 
turbine blades. Fan blades are the most flexible compo- 
nents, especially those in the high bypass ratio commer- 
cial airplane turbofan engines. In case of compressors, 
demands of increasing pressure ratios and weight reduc- 
tion result in highly loaded stages operating in transonic 
inflow. In addition, the front stages are subjected to se- 
vere inlet flow distortions. 

The nature of flutter is complex in turbomachines, 
and the fundamental mechanisms are often difficult to 
distinguish. Conventionally, the various flutter regions 
of a compressor are demarcated as shown in Fig. 24. 
The mass flow is a corrected value, i.e. an equivalent 
mass flow under ambient conditions. The pressure ratio 
is the ratio of exit stagnation pressure to inlet stagna- 
tion pressure of the compressor. The surge line is the 
highest operating line above which the airflow can re- 
verse abruptly creating rotating stall, called surge. The 
subsonic and supersonic stall flutter are two of the more 
distinct regimes. The subsonic stall flutter occurs in part- 
speed conditions when the inlet flow is high subsonic or 
transonic. Typically, they are a Single Degree of Free- 
dom (SDOF) torsion mode flutter. The supersonic stall 
flutter occurs at high-speed operating conditions when 
the inlet flow is supersonic and contains detached lead- 
ing edge shocks. Often, they are a SDOF bending mode 
flutter. A SDOF flutter, unlike a classical frequency coa- 
lescence type bending-torsion flutter involves energy ex- 
change between a single mode and the fluid. The driving 
phenomena in compressors are either oscillating shock 
waves or stalled flow. A simple linear structural model is 
enough as long as the aerodynamic damping is correctly 
evaluated. The complexity mainly lies in the latter task. 


CFD/CSD Predictions 

The 1990s saw the development of large scale paral- 
lel solvers for unsteady aerodynamics in turbomachines, 
see for example Refs. [93, 94, 95]. The high computa- 
tional requirements restricted the simulations to a single 
stage. More recently, solvers with massive parallelism, 
i.e. with demonstrated scalability on 100s of CPUs have 
been developed [96, 97]. For example, the Stanford Uni- 
versity code TFL02000 has been tested for scalability on 
1024 processors for a model problem with 9 whole annu- 
lus blade rows with 93. 8M grid points. There have been 
limited efforts however, to couple it subsequently to tur- 
bomachinery CSD, see Doi [98]. Successful coupling with 
CSD was carried out later at the University of Maryland 
for helicopter main rotor loads and acoustics prediction 
(see section on rotorcraft aeroelasticity) . 

Researchers have tried applying large scale NS CFD, 
i.e. whole-annulus, multi-blade row, unsteady simula- 
tions, to forced response and flutter since the early 2000s. 
These were uncoupled simulations — either predefined 
blade motions were prescribed to CFD to study the flow 
distortions [99, 100, 101], or CFD calculated flow distor- 
tions were imposed on the structure to calculate blade 
deformations [102, 103]. The first procedure can be used 
for an approximate determination of flutter. The second 
procedure can be used for an approximate determina- 
tion of forced response. The first procedure exploits the 
assumption that flutter in turbomachinery is predomi- 
nantly a SDOF phenomena. Thus prescribing the torsion 
mode (or bending mode for supersonic flutter) at its nat- 
ural frequency, and then integrating blade airloads over 
each cycle to calculate damping is accepted as a reason- 
able approach. 

Pioneering work on coupled analysis began during 
the 1990s used frequency domain coupling, see for exam- 
ple Williams et al. [104] for panel solvers, and Geroly- 
mos [105] for Euler solvers. In the work by Geroly- 
mos, the structural dynamics (for a particular mode) was 
solved in the frequency domain. The periodic response 
was then transferred to the 3-D Euler solver. The Eu- 
ler solver was marched in time domain, until periodicity. 
This cycle of periodic CSD and time domain CSD was 
repeated until convergence. Note that this approach, if 
applied to a helicopter rotor, will lead to rapid diver- 
gence. A similar method, but with a key innovation, is 
used in helicopter rotors. This method, termed the delta 
method, is discussed later in the section on rotary wing 
aeroelasticity. 

Time domain coupling using staggered algorithms 
began to be applied in the 2000s with single sector cal- 
culations for a single blade-row. Carsten et al. [106] have 
systematically compared predicted modal decay between 
the time domain and frequency domain solvers in the 
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(a) The first bending and first torsion eigen- 

modes (b) Comparison of logarithmic decrement for frequency and time domain 

method, first bending(left) and first torsion(right) 


Figure 25: RANS/FEM flutter calculations for a 4-bladed segment of a low pressure compressor annulus 
using 5 structural modes (20 modes in total); 8 Pentium-4 CPUs; Carsten et al. 2003 
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(a) Comparison of two fan intake geometries 


(b) Least stable mode vs % speed 


Figure 26: A RANS/FEM study of intake duct effects on fan flutter stability; Vahdati et al. 2002 


presence of oscillating shocks. For the conditions studied 
no significant differences were found, as shown in Fig. 25. 
Whole annulus calculations have been performed exten- 
sively by researchers at the Imperial College, London. 
The original work by Imregun et al. [107, 108] showed 
that forced vibration predictions remained unchanged 
between a typical sector analysis (even with adjustment 
of blade numbers) and whole annulus calculations, as 
long as multi blade row interactions were not impor- 
tant. Subsequent work by Sayma et al. [109], Breard 
et al. [110], and Vahdati et al. [Ill] focused on multi 
blade row simulations. The last study showed the im- 
pact of intake duct geometry on the flutter stability of 
inlet fans, see Fig. 26. Flutter and forced response calcu- 
lations for core compressors are more involved, compared 
to fan flutter, because of the close spacing of blade rows. 


Rotor-stator blade row interactions involve multiple up- 
stream and downstream rows. In 2005, Wu et al [112] 
demonstrated the initial feasibility of a massive RANS 
simulation (Baldwin-Barth turbulence model) including 
17 blade rows, see Fig. 27(a). It was performed on 128 
CPUs on SGI Origin 3000 with a relatively coarse (com- 
pared to the expected requirements) grid of 68M points. 
A more detailed study on a 5 blade row model was car- 
ried out recently by Vahdati et al. [113], see Figs. 27(b) 
and 27(c). The structural model in all of the above stud- 
ies was elementary in comparison to CFD. Only rotors 
1 and 2 were modeled, and only the 1st flap and 1st 
torsion modes (of various nodal diameter families) were 
considered. 

Fleeter and his co-researchers have performed fully 
coupled FEM based CFD / CSD simulations of turboma- 
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0) (c) 

Figure 27: A HPC based Whole-annulus, multi blade row aeroelastic analysis of axial-flow compressors 
using unstructured RANS / FEM coupling; (a) A core-compressor coarse 17 blade row mesh; (b) A 
refined 5 blade row mesh; (c) Structural mode 11, susceptible to 32nd engine order excitation; Wu et 
al. 2005; Vahdeti et al. 2007 


chinery flows both for flutter [114] and for forced re- 
sponse [115]. The latter work involved the vibratory 
stress predictions of a 1.5-stage axial flow compressor 
(Purdue Transonic Compressor) for a inlet guide vanes- 
rotor-stator configuration of 18-19-18 blades. Only the 
rotor and stator were modeled. One segment from each 
blade row was placed in one CPU. The steel stator vanes 
were modeled using 2464 hexahedral brick finite ele- 
ments, 4 across chord and 14 across span. The third 
mode frequency was 5203 Hz, close to the part-speed 
blade passage frequency of 4750 Hz (part-speed operat- 
ing RPM is 15,000, number of blades 19). A simulation 
was performed at 5203 blade passage frequency (with 
RPM 16,430). First, the unsteady periodic flow of the 
rotor was calculated. This required 31 rotor revolutions 
(80h of CPU time on 2 IBM SP2 processors). The stator 
vanes were not allowed to deform during this time. This 
was an aerodynamic only simulation for which a large 
time step (3844 per period) could be taken. Next, the 
fully coupled transient fluid-structure simulation was ini- 
tiated. For this case, the time step was reduced by three 
times. The solution was marched until steady state pe- 
riodic conditions. This required 700h of CPU time on 2 
IBM-SP2 processors. 

ROTARY WING AIRCRAFT 

A recent review of rotary wing CFD/CSD, detailing 
the historical evolution and current status of airloads and 


structural loads prediction as of 2006, can be found in 
Datta et al. [116]. The following sections present a sum- 
mary of that review. In addition, the reader is brought 
up to date with the latest developments. 

As before, the focus here is on CSD, and coupled 
high-fidelity CFD/CSD methods. A thorough review of 
rotorcraft CFD methods, from its early inception to its 
evolution into its present status to future challenges, can 
be found in Strawn et al. [117] in 2007. As before, the 
review is divided into two sections. The first section sum- 
marizes the status of detailed CSD. The second section 
summarizes the key aeroelastic phenomena, and the sta- 
tus of CFD/CSD in predicting them. 

Structural Modeling 

The non-linear coupled PDEs of rotating blade 
dynamics were systematically derived in the 1970s 
by several researchers, notably Hodges and his co- 
researchers [118, 119], Kvaternik et al. [120], and Rosen 
and Friedmann [121]. The rotor blade was treated as 
isotropic and straight Euler-Bernoulli beams. By retain- 
ing the second-order effects of elastic motion in the ki- 
netic and strain energy terms, and by adding the non- 
linearities (up to second order) in extension and torsion 
produced by bending, ‘almost-exact’ beam models can be 
developed. For larger bending deformations, the nonlin- 
earities can be handled by using the ‘geometrically exact’ 
beam theories [122], or within a multibody type formu- 
lation [123] which allow both rigid and elastic motion of 
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structural components. 

Since the mid 1980s, there was an emphasis on im- 
proving cross-section analyses to account for anisotropy 
produced by composite materials. Research in compos- 
ite blade modeling has been reviewed in detail by several 
researchers, most recently by Friedmann in 2004 [124]. 
The basic idea has been to perform a linear or non-linear 
two-dimensional cross-sectional analysis, uncoupled from 
the ID non-linear beam theories mentioned above, with 
the goal of obtaining the cross-sectional elastic constants 
and warping. It has been shown using asymptotic meth- 
ods that such uncoupling assumptions hold true in most 
cases for rotor blades [125]. The cross-sectional elastic 
constants defining the stress-strain relation are however 
fully coupled. For low frequency dynamic behavior, a 
classical analysis with a single shear strain measure (as 
in the case of isotropic beams) can be reasonably accu- 
rate, as long as the coupled elastic constants are calcu- 
lated properly [126]. For higher accuracy, as required for 
rotor analysis, the strain vector can include two addi- 
tional strains to account for the transverse shears, or one 
additional strain to account for the Vlasov or restrained 
warping. The accuracy of the model depends more on 
the accuracy of the elastic constants than the number of 
strain measures. There are several methods for extract- 
ing these constants. They are: engineering methods that 
work well for certain type of cross-sections [127], asymp- 
totic methods which yield closed form solutions for thin 
walled beams [128, 129], and finite element methods. The 
latter can be based either on St. Venant’s principle [130] 
or on asymptotic methods [131]. The latter reference pro- 
vides a generalized treatment, where the cross-sectional 
equations and the one dimensional blade equations are 
rigorously deduced from three-dimensional elasticity the- 
ory. These analyses or static measurements are used to- 
day to obtain the section properties required to operate a 
beam-theory based comprehensive code. The dynamics 
of the blade is then calculated, subject to proper model- 
ing of the boundary conditions. The latter is a challeng- 
ing task because of nonlinearities associated with the root 
end of the blades, the prime examples of which are lag 
dampers and control systems. 

Gandhi and Chopra [132] incorporated an analytical 
elastomeric lag damper model in a rotor comprehensive 
analysis. Finite element based damper models have also 
been developed, such as by Smith et al [133]. In general, 
the difficulties with modeling lag dampers are numerous. 
Non-linearities associated with temperature, amplitude, 
and frequency of these structures create unknown param- 
eters and reduce prediction accuracy of the blade models. 
Complex boundary conditions can be handled in a gener- 
alized manner using a multibody formulation [134]. Orig- 
inally developed for spacecraft dynamics, it provides a 
framework to treat rigid and flexible structural elements 


undergoing arbitrary rotations and translations with re- 
spect to one another. In addition, unlike conventional 
formulations, multibody formulations enable the increase 
in scope of analysis without having to re-derive and re- 
validate the equations with each additional feature. One 
of the early rotor aeromechanical stability analysis devel- 
oped to represent arbitrary rigid and flexible dynamics 
was GRASP by Hodges et al. [135]. This multibody like 
analysis was designed to treat modern rotors with ar- 
bitrary boundary conditions in steady, axial flight and 
ground contact conditions. Multibody formulations en- 
able large deformation problems to be accurately treated 
while still retaining second-order nonlinear beam theory 
(i.e., without using geometrically exact theories). This is 
done by simply breaking up the rotor blade into multi- 
ple bodies each undergoing only moderate deformations 
in its local frame. The total deformation is obtained by 
adding the local deformations to the end contribution of 
the adjoining parent element as a rigid motion. 

DYMORE [136] and MBDyn [137] are examples of 
modern multibody structural codes. These codes are of- 
ten geared to solve the full non-linear finite element equa- 
tions using time-marching schemes. This is due to the 
presence of algebraic constraint equations. Time march- 
ing schemes are less suitable for rotor trim. Artificial 
damping is often required to suppress the natural re- 
sponse of the lowly damped modes, particularly the lag 
mode. The damping needs to be subsequently removed 
as the solution progresses to arrive at the correct steady 
state periodic response. In addition, the use of algebraic 
equations prevents the use of state space based aeroelas- 
tic stability analysis. The accuracy of blade loads predic- 
tion, however, have not shown any significant improve- 
ment, so far, with the expansion of the scope of structural 
modeling. The importance of detailed boundary condi- 
tion modeling, and their impact on the current deficien- 
cies in blade loads prediction, remains to be understood. 

The development of coupled rotor-fuselage dynamic 
analysis can be traced back to the early work of Gersten- 
berger and Wood [138] in 1963. The work was based on 
a method of impedance matching, which was later used 
extensively during the 1970s and 80s [139, 140, 141, 142]. 
Beginning 1984, NASA Langley Research Center carried 
out a successful Design Analysis Methods for Vibrations 
(DAMVIBS) program, with the active participation of 
four major rotorcraft manufacturers (Bell, Boeing, the 
former McDonnell-Douglas, and Sikorsky). As part of 
the program detailed FEM of several metal and compos- 
ite helicopter fuselages were constructed, e.g. the AH- 
1G, Model 360, AH-64A, OH-6A, CH-47D, S75, D292, 
and the UH-60A. An overview of the program can be 
found in Kvaternik [143] . The description of the airframe 
FEM models as developed by each of the companies can 
be found in Refs. [144, 145, 146, 147]. The prevalent 
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methods of rotor-fuselage dynamic analysis were sum- 
marized by Kvaternik et al. [148]. 

During the 1990s, iterative or fully coupled formula- 
tions began to be applied to the rotor-fuselage problem. 
An early work was that of Stephens and Peters [149], fol- 
lowed by Vellaichamy and Chopra [150], and Chiu and 
Friedmann [151]. During the 2000s, detailed 3-D fuse- 
lage FEM models began to be coupled to rotor dynam- 
ics using these formulations. Cribbs et al. [152] studied 
active vibration control techniques in the fixed frame us- 
ing a MBB BO-105 type four-bladed hingeless rotor he- 
licopter. The DAMVIBS AH-1G fuselage was used by 
Yeo and Chopra [153] to carry out vibration prediction 
studies (see Figs. 6 and 7). A similar (to DAMVIBS) 
SH-60 fuselage was later used by Yang et al. [154], to 
study the vibration impact of damaged rotors. Bauchau 
et al. [155] applied the method of component mode syn- 
thesis to couple fuselage dynamics to a multibody rotor 
model. 

Except for the two-bladed teetering AH-1G heli- 
copter, the inclusion of fuselage dynamics had a minor 
effect on the blade loads. Predicted vibration in the fuse- 
lage, of course, relies on fuselage dynamics. As mentioned 
earlier, all of these studies focused on vibratory forcing 
transmitted via the shaft. The effect of rotor wake im- 
pingement were not included. 

Aeroelastic Phenomena and CFD/CSD Predic- 
tions 

In rotary wing aeroelasticity, a first principles solu- 
tion requires more than CFD and CSD. It requires, in 
addition, a method for simultaneously predicting the un- 
steady control pitch angles and aircraft operating states. 
This is referred to as the Vehicle Flight Dynamics (VFD) 
problem in this paper. The requirement arises be- 
cause unlike fixed-wings and turbomachines, the dynamic 
loads, flutter (flap-lag), and air-resonance characteristics 
of the rotor blades are coupled inherently and nonlin- 
ear ly to vehicle dynamics via the control pitch angles 
prescribed at the blade root. 

In steady flight, straight-level or turns, VFD deter- 
mines the operating state which holds the entire aircraft 
in equilibrium. The operating state is defined by the 
rotor/s pitch control angles (only collective in the case 
of tail rotors), aircraft attitude angles and rates, fuselage 
properties (aerodynamic interaction, weight-distribution, 
tail properties and setting, compounding and special con- 
figurations). The word ‘coupled trim’ is used by dynam- 
icists to describe the VFD associated with steady flight. 
During unsteady flight, VFD is a generalization of cou- 
pled trim, and determines the operating state which fol- 
lows an intended flight dynamics of the aircraft. For a 
general maneuver, these intended targets can be a pre- 


scribed trajectory in space, or a prescribed mean loading 
schedule (e.g. g-level and a prescribed set of hub roll 
and pitch moment variation with time). The optimal 
control problem of determining the control angles based 
on these intended targets is still beyond state-of-the-art, 
even with lower order free wake based lifting-line aero- 
dynamic models. 

For steady flight, the VFD problem is relatively sim- 
ple. The trajectory is trivial, and the mean loading 
schedule can be calculated from first principles using a six 
DOF equilibrium in level flight, or an eight DOF equilib- 
rium in a most general steady maneuver (a coordinated 
helical turn). Calculating VFD (or ‘coupled trim’) is the 
first step for predicting performance, loads, vibration, 
stability, and acoustics. An innovative method for si- 
multaneous determination of VFD during CFD/CSD for 
rotors was described by Tung et al. in 1986 [156]. Re- 
ferred to in this paper as the ‘delta coupling’, it has re- 
mained the most efficient approach and has been the cor- 
nerstone of all current advances in rotorcraft CFD/CSD. 
The method is also referred to as ‘loose coupling’, how- 
ever, it bears no resemblance with CFD/CSD loose cou- 
pling in fixed wings, and the confusion is best avoided. 

For unsteady maneuvers, the analysis becomes con- 
ceptually simpler when the rotor control angles and ve- 
hicle dynamics are known. The VFD is then uncoupled 
from CFD/CSD simulation. Prescribing VFD, either 
from measured data or from isolated lower order pre- 
dictions, reduces the CFD/CSD problem to a straight- 
forward, partitioned, numerical integration problem, 
similar to those applied in fixed wing and turbo machin- 
ery disciplines. Except that, in rotorcraft, the problem 
is currently solved at a much lower fidelity compared to 
the other disciplines. First, the structural model for an 
isolated rotor blade is almost always a beam, compared 
to 3-D structures in other disciplines. Second, the em- 
phasis in rotors has been on blade loads, not flutter, and 
therefore errors in approximate fluid-structure coupling 
methods and low-order (below second order) temporal in- 
tegration have been of lesser consequence. There is only 
a single example of CFD / CSD coupling performed so far 
(see later) for a realistic unsteady maneuver. 

The application of delta coupling during the 1990s 
yielded limited payoffs compared to lifting-line meth- 
ods [157, 158, 159]. A comprehensive report on loads 
prediction for the Research Puma helicopter using U.S. 
and European CFD/CSD methods of the time is docu- 
mented in Bousman et al. [160] in 1996. Philippe [161], 
in 1994, presented a summary of other European efforts. 
The key drawback at the time was the lack of consis- 
tent coupling due to the inability of the CFD solvers in 
resolving the grid and compressibility effects. Renewed 
interest in the 2000s was spurred by the emergence of 
reliable unsteady Euler/RANS solvers. 
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The resurgence of rotary wing CFD/CSD/VFD, in 
the U.S., stems from the systematic work of Datta et 
al. [162, 163] and Potsdam et al. [164]. The advances were 
enabled in substantial part by the reliable and repeat- 
able flight test data provided by the NAS A/ Army UH- 
60A Airloads Flight Test Program. Reference [162] dis- 
sected the coupled aeroelastic loads problem into aerody- 
namic, structural dynamic, and trim components; identi- 
fied the key mechanisms of rotor vibration at high speed, 
isolated the effects of wake and transonic pitching mo- 
ments, and clarified the key contribution of Euler/RANS 
over lifting-line models. Reference [163] brought the 
components back together again and demonstrated im- 
proved airloads and structural loads prediction using 
CFD/Comprehensive Analysis coupling. The predictions 
were compared systematically between lifting-line air- 
loads, CFD airloads, and measured airloads. The CFD 
airloads were computed using near blade RANS coupled 
to free-wake inflow using the field velocity approach by 
Sitaraman and Baeder [165]. 

Potsdam et al. [164] focused on airloads. However, 
the work extended RANS capabilities beyond the high 
speed flight to a wake dominated transition flight and a 
stall dominated highly loaded flight. In the wake dom- 
inated flight, the work demonstrated that CFD is capa- 
ble of predicting the dominant wake induced vibratory 
airloads at transition speeds even without fully captur- 
ing the individual vortices in the wake. At the dynamic 
stall flight, the work demonstrated that a one equa- 
tion turbulence model (Baldwin-Barth) is adequate for 
predicting the basic stall phenomena on the retreating 
side. The stall flight was subsequently studied in de- 
tail by Datta and Chopra [166], isolating the physics of 
structural dynamics, wake, stall, and trim. The physics 
of the stall mechanisms and their impact on predicted 
pitch-link loads using CFD / Comprehensive Analysis cou- 
pling were clarified. The CFD analysis used in refer- 
ences [162, 163, 166] was subsequently extended to in- 
clude overset capabilities for wake capture, similar to 
that of reference [164], by Duraisamy and Baeder [167]. 
This methodology was then used to carry out predictions 
of vibratory airloads and structural loads at all three level 
flight conditions in Refs. [168, 169]. The work of Pots- 
dam et al. [164] was subsequently extended by Lim et 
al. [170] for studying vortex loadings in more details. The 
focus was on descent flights of the UH-60A and HART 
II rotor where, both the vortex induced loadings as well 
as the individual BVIs are present in the airloads. The 
BVI airloads predicted by CFD were improved compared 
to lifting-line predictions. A significant portion of the 
improvement was in low frequency harmonics (3/rev), 
brought about by vortex loadings in the first and fourth 
quadrants, rather than individual BVI. Increasing the 
background grid resolution from 10% chord to 5% chord, 


a difference of 19. 4M and 107. 6M grid points, did not 
show any improvement in the individual BVI loadings. 

Concurrent developments in Europe over the last 
five years have also contributed to the increased un- 
derstanding of CFD/Comprehensive Analysis capabili- 
ties for loads prediction in helicopter rotors. The high 
speed wind tunnel test airloads of the ONERA 7A and 
7AD rotors have been studied and predicted by sev- 
eral researchers, see Servera et al. [171], Altmikus et 
al [172], Pomin and Wagner [173], and Pahlke and Van 
Der Wall [174]. A discussion on the key accomplishe- 
ments in all of the above work can be found in the review 
article, Ref. [116]. 

The developments cited above provided a renewed 
impetus for integrating CFD reliably within the rotor 
design process. At the time, among all key design is- 
sues, noise, tied directly to vulnerability, was considered 
to be of immediate urgency. Innovative blade shapes 
and advanced blade geometries for noise reduction are 
not considered as serious contenders for blade design be- 
cause their impact on performance, loads, and vibration 
cannot be predicted with confidence. As CFD begins to 
predict the fundamental mechanisms for structural loads 
— 3-D unsteady transonic pitching moments, wake in- 
duced airloads and BVI, and dynamic stall cycles — a 
considerable attention is now focused on increasing CFD 
throughput while maintaining the same level of accu- 
racy as benchmarked by the above cited works. The 
DARPA Helicopter Quieting Program was designed to 
carry out such a task. A key objective of the DARPA 
program was to adapt non-rotorcraft CFD with demon- 
strated scalability to rotor airloads prediction. For exam- 
ple, the TFLOW2000 code was adapted for rotary wing 
CFD (extended and renamed SUmb) and coupled to a 
rotorcraft comprehensive analysis, as part of the pro- 
gram. An additional objective of the program was to 
investigate the impact of high fidelity wake and turbu- 
lence modeling, e.g. methods for solution adaptive back- 
ground mesh refinements, different near-blade and back- 
ground solvers, and v 2 -f and KES turbulence models, 
which these analyses incorporated. The focus was not on 
CSD or coupling methods. To this end, trimmed aeroe- 
lastic blade deformations from Ref. [164], were supplied 
to the participants as an initial step. As a follow on, the 
participants performed the same CFD/CSD/VFD cou- 
pling as in Refs. [163, 164] using the same comprehen- 
sive analysis tools. The validation set was extended to 
include acoustic signatures in addition to performance 
and airloads. The DNW model UH-60A, and the BO- 
105 model HART-II rotor were considered as valida- 
tion cases for acoustics. The contributions of the three 
program participant teams are partially documented in 
Refs. [176, 177, 178, 179], 

All of the above coupled studies, with the exception 
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of reference [173], used variations of the delta coupling 
method in the time domain. Reference [173] does not 
calculate control angles as part of the CFD coupled so- 
lution but determines them a priori using an uncoupled 
lifting-line comprehensive analysis. The delta method is 
unique to rotorcraft and is not found in fixed-wing or 
turbomachinary aeroelasticity. It accounts for the fol- 
lowing key requirements of rotor CSD: (1) it operates, 
by design, at or near resonance in flap, (2) it requires 
simultaneous solution of VFD for the pitch control an- 
gles. The description of the original method can be 
found in Ref. [156], its current adaptations in any of 
the above references, see for example Refs. [164, 163]. 
The method is implied in rotorcraft literature by the fre- 
quent use of the term ‘CFD/ Comprehensive Analysis’. 
The term ‘Comprehensive Analysis’, by itself, histori- 
cally refers to a complete analysis, which contains within 
itself all of the airloads, structural loads, and VFD mod- 
els. Thus, ‘CFD/ Comprehensive Analysis’ simply refers 
to a CFD/CSD/VFD analysis which uses the CSD and 
VFD of existing comprehensive analysis while replacing 
the built-in lifting- line/surface aerodynamics with CFD 
via the delta coupling method. 

In the special case where the control angle varia- 
tion is assumed to be known a priori, either from mea- 
sured test data or from separate uncoupled estimates 
from lower order models, the problem reduces to a pure 
CFD/CSD analysis. A straight-forward, partitioned, nu- 
merical integration procedure can then be applied, sim- 
ilar to those used in fixed-wing or turbomachinery. Ex- 
cept that, for an isolated rotor analysis, the problem is 
rendered much simpler as the structural model is almost 
always a beam, and a low order fluid-structure interface 
is much simpler to devise. These methods are known was 
‘tight coupling’ in rotorcraft. Again, note the confusing 
nomenclature as compared to fixed wing terminology. In 
fixed wing tight coupling implies strict time accuracy via 
sub- iterations. One of the first implementations of tight 
coupling in rotorcraft using RANS CFD was by Bauchau 
and Ahmad [180]. Recent studies are those of Pomin and 
Wagner [173] and Bhagwat et al [181]. The latter work 
is of significance as it is the first application of CFD to 
a real unsteady maneuver. The method is applied here 
to a pull-up maneuver, the second most severe (loads) 
of the UH-60A maneuvers. Note that, the advances in 
this study were enabled by the NAS A/ Army UH-60A 
Airloads Flight Test Program Data. The rotor control 
angles and vehicle dynamics were prescribed from mea- 
sured flight data. Several key conclusions can be drawn 
from this work. First, it demonstrates RANS capability 
for predicting advancing blade transonic stall. Second, 
it appears to negate the role of fuselage upwash on ro- 
tor stall and pitch-link loads during a pull-up maneuver. 
Third, it demonstrates that an isolated rotor calculation 


can predict the oscillatory (peak-to-peak) control loads 
correctly in the severest of the UH-60A maneuvers. 

A LIST OF TEMPORAL COUPLING 
PROCEDURES 

In view of the conflicting nomenclatures often used 
across disciplines to describe the various methods of tem- 
poral coupling, a list is summarized below. 

1. Loose coupling: Partitioned formulation based tem- 
poral advancement, but without sub-iterations. 
Predictor-corrector schemes often used to seek near- 
second order accuracy. Standard nomenclature in 
fixed wing. Described as tight coupling in rotary 
wing. 

2. Tight coupling: Partitioned formulation based time 
advancement, but with sub-iterations. Temporal ac- 
curacy level of worse solver. Standard nomenclature 
in fixed wing. Also called close coupling or strong 
coupling. 

3. Periodic coupling: For periodic dynamics. First ap- 
plied by Gerolymos for turbomachinery. Iterates on 
harmonics of loads and deformations. Strict time 
accuracy can be enforced on any chosen set of har- 
monics. 

4. Delta coupling: Periodic coupling with the use of 
aerodynamic damping from lower order models. 
Only successful approach for CFD/CSD/VFD in ro- 
torcraft. Conflictingly nomenclature as loose cou- 
pling in rotary wing. Note, not loose in the fixed 
wing sense of the term. Strict time accuracy can be 
enforced on any chosen set of harmonics. 

5. Full coupling: Formulation coupled via boundary 
conditions (deflection and velocity compatibility, 
and stress equilibrium) and integrated together. 
Also called direct coupling , intimate coupling, or 
monolithic coupling. 

CONCLUDING OBSERVATIONS 

The objective of the present paper was to moti- 
vate research towards the next generation HPC based 
high-fidelity rotorcraft comprehensive analysis. To this 
end, the frontiers of fixed wing, turbomachinery, and ro- 
tary wing aeroelasticity were reviewed, along with a sur- 
vey of the fluid-structure coupling methodologies used 
in each. The underlying theme of the paper was rotor- 
craft. Therefore, the technology areas in each discipline 
were tied to the future rotorcraft requirements, where 
deemed relevant by the authors. In addition, compar- 
isons were drawn between the computational procedures 
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used in each of the disciplines, and those required by the 
unique rotorcraft requirements. 

CFD is only one part of a next generation com- 
prehensive analysis. The ultimate objective of reliably 
discriminating between innovative blade shapes, rotor 
types, and rotor-fuselage configurations, requires CFD, 
CSD, high-fidelity coupling, and a solution procedure 
that ties CFD/CSD to VFD satisfying the unique oper- 
ating principles of rotorcraft. The focus so far has largely 
been on CFD. This is because CFD, up to recently, was 
unable to meet the fundamental requirements for vibra- 
tory structural loads in rotorcraft: predicting transonic 
pitching moments on blades, wake induced airloads on 
both blades and fuselage, dynamic stall cycles on blades, 
and buffet loads on the fuselage driven by massive sep- 
aration. HPC based CFD is now beginning to address 
these requirements. The methods of high-fidelity CSD, 
fluid-structure coupling, and solution procedures must 
now be re-evaluated and re-shaped for taking rotorcraft 
comprehensive analysis to the next level. 
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